雑に読む「Numerical Recipes in C 日本語版」
https://m.media-amazon.com/images/I/512QW6W25KL._SX355_BO1,204,203,200_.jpg#.png https://www.amazon.co.jp/dp/4874085601
1993
動機
2022後期に数値計算のレポートがあったのだが、数値計算になれてなさすぎたせいか、何もアイデアが浮かばなかった 数値計算の手法をざっくり把握しておきたい
目標takker.icon
継続
何日かに1回読む
井戸端を開いたときに必ず書き込む
ちょっと強すぎるかも？takker.icon
分厚い本なので、机に座らないと読めない
table:読んだ日
2023-03-20
2023-03-22
コードはJSで書く
Python版も書きたいyosider.icon
Python苦手～><takker.icon
通読だけなら、7章までざっと読んである
具体的に数式やコードを書いて動かしたわけではないtakker.icon
少なくとも、2章は数式やコードまで全部追いたい
完了
一通りコードを書いて動かしたいが、分量が多すぎて終わりそうにない
数式の理解は飛ばしていい
分量が多すぎて終わらない
第2章あたりの計算はさっとコードを書けるようになる
数値計算コードを書くのに慣れたい
関数の計算あたりは考え方がわかればいい
いつかやる
Rustで書き直す
第１章 準備
1.0 はじめに
1.1 プログラムの構成と制御構造
知ってるので略takker.icon
(説明できるとは言っていない)
1.2 科学計算のためのC言語プログラム作法
C言語の地雷のことしか書かれてないので略takker.icon
1.3 誤差，正確さ，安定性
2023-02-16
ここから書いてくtakker.icon
intやlongは固定小数点としている
整数$ \subseteq固定小数点ということだろうかtakker.icon
まあ間違ってはいない
$ \mathsf{value}=s\cdot M\cdot B^{e-E}
ここでBは基数(基底)で,通常B=2であるが,B=16のこともある
$ B=16の処理系なんてあるのかtakker.icon
2023-02-26
$ s=\pm1
$ 0\le M<2^{23}=8388608
32bit符号付き整数の最大値が$ 2^{31}だから、正確に表せる数値の範囲が整数より狭いだいぶ狭いことがわかる
$ 0\le e-E<2^8=256
10進数に換算すると何桁くらい？takker.icon
$ 2^m={10}^n
$ \iff m\ln2=n\ln10
$ \iff n=m\frac{\ln2}{\ln10}=0.30102999566\cdots\times m
$ \therefore B^{e-E}<2^{256}\simeq{10}^{77.06}
77桁くらいか
あれ？256のうちどのくらいが負の指数に使われるんだ？
指数の下駄$ Eがそれを決めているのかな？
table:floatの例
s 指数 仮数
1/4 0 01111111 1000000000000
1/2 0 10000000 1000000000000
3 0 10000000 1100000000000
はい？なんで$ 1/2=1\times2^{-1}と$ 3=3\times2^0の指数部の値が同じなの？？takker.icon
$ 1/4=1\times2^{-2}の指数部が1/2の指数部を反転した値になっているのも謎
浮動小数点表現では, ビットの並びが異なっても同じ数値を表すことがありうる。
例えば B = 2 で, 仮数部 (浮動小数点表示の) の左側 (高位の桁) 0のビットがある場合は、仮数部全体を左にシフトして(つまり2の累乗を乗算して)、その分だけ指数部 (浮動小数点表示の)を減らしても,同じ数値を表す。
この操作により仮数部の最も左側に0のビットが来ないようにすることを正規化 (normalization) という｡
このあたり、数値計算のアルゴリズムとはあまり関係ないから、飛ばして2章を先にやったほうがいい気がしてきたtakker.icon
2023-03-17 13:11:50 飛ばします
いずれ読む
第２章 連立１次方程式の解法
ここから開始
2.0 はじめに
特異な連立方程式、非特異な連立方程式
連立方程式でうまく解が求まらない場合がある
理論的に求まらない場合
退化している方程式は特異(singular)であるという。
連立方程式の場合は、以下に分けられる
1つ以上の方程式が、他の方程式の線型結合で表されるもの
例
$ \begin{pmatrix}1&2&3\\7&-1&4\\8&1&7\end{pmatrix}\begin{pmatrix}x\\y\\z\end{pmatrix}=\begin{pmatrix}-4\\5\\1\end{pmatrix}
1行目と2行目を足すと3行目になる
$ \begin{pmatrix}1&2&3\\7&-1&4\end{pmatrix}\begin{pmatrix}x\\y\\z\end{pmatrix}=\begin{pmatrix}-4\\5\end{pmatrix}と同値
また、どの方程式にも、いくつかの変数が全く同じ線形結合となって現れる(この条件を列退化 column degeneracyという)ならば、やはり解は一意には定まらない
言っていることがわからないtakker.icon
いくつかの変数が全く同じ線形結合となって現れる
「なんの」線形結合？？
例$ \begin{pmatrix}1&2\\7&-1\\8&2\end{pmatrix}\begin{pmatrix}x\\y\end{pmatrix}=\begin{pmatrix}-4\\5\\1\end{pmatrix}
計算誤差で求まらない場合
いくつかの方程式が非常に線型従属に近い場合
未知数が非常に多いときに特に注意
連立方程式を解くpackageは、これらの病理の検出・防止のせいで複雑になってしまう
大雑把に言って、以下の場合なら単純に計算して問題ないらしい
20~50個の未知数をfloatで解く
数百個の未知数をdoubleで解く
行列
記法の解説
計算線形代数の仕事
2章で扱う問題の紹介
$ A_{ij}x_j=b_i\quad,i\le M,j\le Nのとき、$ M>Nな方程式を優決定方程式と呼ぶらしい これは解が不能にならない場合もそうなる？
例：$ M-N個の方程式が、それ以外の方程式の線型結合で表せるとき
この際厳密解は存在しないが、できるだけ正確に満たす値を最小二乗法で求めることはできる $ \pmb{A}^\top\cdot\pmb{A}\cdot\pmb{x}=\pmb{A}^\top\cdot\pmb{b}
標準的なサブルーチンパッケージ
C言語の場合
linpack以外は化石ライブラリと化していそうtakker.icon
Rustは前調べてたことがある。まだデファクトスタンダードが決まっていなかったはず
JSはあるのかな
スクリプト言語で行列計算する意味を問うてはならない
この辺からprogramの出番が出る
メリット
$ \pmb{A}\cdot\pmb{C}=\pmb{B}から逆行列$ \pmb{A}^{-1}と方程式の解$ \pmb{C}を同時に求められる
algorithmが単純で理解しやすい
デメリット
すべての右辺を同時に格納・操作しなければならない
意味がわからなかったtakker.icon
何と同時？
どこに格納？
$ \pmb{A}\cdot\pmb{C}=\pmb{B}から解$ \pmb{C}を求めるだけなら、これより3倍速い方法があるらしい
連立方程式を解くのにも逆行列を求めるにも、これよりLU分解のほうがいい algorithmが単純で理解しやすいから
列拡大行列の消去法
任意個の連立方程式と逆行列を一つの連立方程式に合体し、一度に解く
任意個の連立方程式$ \pmb{A}\cdot\pmb{x}_i=\pmb{b}_iと$ \pmb{A}\cdot\pmb{Y}=\pmb{I}を組み合わせて、
$ \pmb{A}\cdot(\bigsqcup_i\pmb{x}_i)\sqcup\pmb{Y}=(\bigsqcup_i\pmb{b}_i)\sqcup\pmb{I}
e.g. $ \begin{pmatrix}0\\3\\4\end{pmatrix}\sqcup\begin{pmatrix}5\\7\\-1\end{pmatrix}=\begin{pmatrix}0&5\\3&7\\4&-1\end{pmatrix}
を作る
これを解くと
$ (\bigsqcup_i\pmb{x}_i)\sqcup\pmb{Y}=\pmb{A}^{-1}\cdot(\bigsqcup_i\pmb{b}_i)\sqcup\pmb{I}
$ =(\bigsqcup_i\pmb{A}^{-1}\cdot\pmb{b}_i)\sqcup\pmb{A}^{-1}
となる。
しらない言葉だ
ダンスやバスケットボールにおける片足旋回のこと
関係なさそうtakker.icon
まだ理解してない
CコードをTSに書き直す
うわ要素番号1始まりかよtakker.icon
まあ直すの大して難しくないけど
2次元配列と行列との対応メモ
https://kakeru.app/a4539cfd60a8e915fa27a79c67284218 https://i.kakeru.app/a4539cfd60a8e915fa27a79c67284218.svg
$ bはこの関係が逆転しているので注意
code:mod.ts
export type Matrix = number[][];
export type Vector = number[];
/** 完全pivot選択のGauss-Jordan法で複数の連立方程式の解と逆行列を求める
*
* @param a 連立方程式の係数
* @param bs 求めたい連立方程式の右辺
* @return 連立方程式の解とaの逆行列
export const pivotGaussJordan = (a: Matrix, ...bs: Vector[]): [Matrix, ...Vector[]] => {
// 次元チェック
const n = a.length;
if (!a.every((row) => row.length === n)) throw Error(a must be a square matrix);
const m = bs.length;
if (!b.every((column) => column.length === n)) throw Error(b must be the same column length as a);
/** pivot選択記録用 */
const indxc: number[] = [];
/** pivot選択記録用 */
const indxr: number[] = [];
/** pivot選択記録用 */
const ipiv: number[] = new Array(n).fill(0);
let big = 0;
for (let i = 0; i < n; i++) {
big = 0;
for (let j = 0; j < n; j++) {
if (ipivj === 1) continue; for (let k = 0; k < n; k++) {
}
ここまで読んでる
疑問点
LU分解から秒で求められる
LU分解から秒で求められる
非正方行列でもできる分解
2.11 逆行列の計算は$ N^3乗の過程か？
おまけ節
$ O(N^3)より計算量を減らすことができるが、計算が複雑になるしそこまでする必要ある？というお話takker.icon
第３章 補間と補外
3.0 はじめに
3.6 ２次元以上の補間
第４章 関数の積分
4.0 はじめに
4.2 初等的なアルゴリズム
4.6 多次元積分
第５章 関数の計算
5.0 はじめに
5.1 級数と収束性
5.3 多項式と有理関数
5.5 ２次方程式，３次方程式
5.7 Chebyshev 近似した関数の微分と積分
6.0 はじめに
7.0 はじめに
7.5 データ暗号化に基づく乱数列
8.0 はじめに
8.3 索引づけと順位づけ
第９章 非線形方程式と非線形連立方程式の解法
9.0 はじめに
9.4 微分を利用した Newton-Raphson 法
9.5 多項式の根
第１０章 関数の最大・最小
10.0 はじめに
10.3 １階導関数を使う１次元探索
11.0 はじめに
11.3 ３重対角行列の固有値，固有ベクトル
11.7 逆反復法による固有値の改良と固有ベクトルの計算 第１２章 フーリエ変換
ここまでやると、来年度の講義の予習になるtakker.icon
12.0 はじめに
12.6 FFTによる最適(Wiener)フィルタ
12.7 FFTによるパワースペクトルの推定
12.9 時間領域でのディジタルフィルタ
12.11 ２次元以上の FFT
第１３章 データの統計的記述
13.0 はじめに
13.2 メディアンの効率的な探索
13.3 連続データの最頻値の推定
13.4 二つの分布が同じ平均値・分散を持つかどうかの検定
13.5 二つの分布は異なるのか？
13.7 線形相関
13.8 ノンパラメトリック(順位)相関
14.0 はじめに
14.2 データの直線への当てはめ
14.3 一般の線形最小2乗法
14.4 非線形モデル
第１５章 常微分方程式の数値解法
15.0 はじめに
16.0 はじめに
16.6 内点での境界条件，特異点の処理
第１７章 偏微分方程式
17.0 はじめに
付録A 参考文献
付録B プログラムの依存関係
付録C プロトタイプ宣言一覧
付録D ユーティリティルーチン
付録E 複素数演算パッケージ
索引
掲載プログラム販売のお知らせ
訳者あとがき
観客席
辞書を通読するのは自分には難しいのでどう読むのか気になる基素.icon
これ辞書だったんですねtakker.icon
そういえばタイトルにレシピと書いてある